%ssget should be installed.
index = ssget;
f=[224,22,159,1236];
err =[sqrt(eps),10^(-5),10^(-14),sqrt(eps)];
rows = index.nrows (f);
index.Name(f)

T0 = 0;
TF = [10,2,2,2];
p = 1;
NumExp = 100;
cputimes = zeros(1,NumExp);
phiptime = zeros(size(f));
phiperr = zeros(size(f));
uphiptime = zeros(size(f));
uphiperr = zeros(size(f));
phipmtime = zeros(size(f));
phipmerr = zeros(size(f));
phipm10time = zeros(size(f));
phipm10err = zeros(size(f));
umphipmtime = zeros(size(f));
umphipmerr = zeros(size(f));
uphipm1time = zeros(size(f));
uphipm1err = zeros(size(f));
uphipmtime = zeros(size(f));
uphipmerr = zeros(size(f));
statephip = zeros(length(f),4);
stateuphip = zeros(length(f),4);
tkchoicestate = zeros(length(f),4);
tkchoicestate10 = zeros(length(f),4);
stateum = zeros(length(f),4);
tkchoicestate1 = zeros(length(f),4);
tkchoicestateu = zeros(length(f),4);


k = 0;

for ff = f
    k = k+1
    Problem = ssget (ff,index);
    sym = index.numerical_symmetry(ff) == 1;
    N = rows(k);
    K = spdiags(ones(p, 1), 1, p, p);
    A = [Problem.A, ones(N, p); zeros(p, N), K];
    ep = zeros(p,1); ep(end) = 1;
    b = [ones(N, 1); ep];
    [wexact, tvals, mv] = expmv_tspan(A, b, T0, TF(k), 100, 'double');
    wexact = wexact(1:N,end);
    
    for i = 1:NumExp
        tic;
        [wp, stats] = phip(TF(k),Problem.A, ones(N, p+1),err(k), sym);
        cpu1 = toc;
        cputimes(i) = cpu1;
    end
    statephip(k,:) = stats;
    phiptime(k) = mean(cputimes);
    phiperr(k) = norm((wexact-wp)./wexact);
    
    for i = 1:NumExp
        tic;
        [wp, stats] = uphip(TF(k), Problem.A, ones(N, p+1),err(k), sym);
        cpu1 = toc;
        cputimes(i) = cpu1;
    end
    stateuphip(k,:) = stats;
    uphiptime(k) = mean(cputimes);
    uphiperr(k) = norm((wexact-wp)./wexact);
    
    for i = 1:NumExp
        tic;
        [wp, stats] = phipmactual(TF(k), Problem.A, ones(N, p+1),  err(k), sym);
        cpu1 = toc;
        cputimes(i) = cpu1;
    end
    tkchoicestate(k,:) = stats;
    phipmtime(k) = mean(cputimes);
    phipmerr(k) = norm((wexact-wp)./wexact);
    
    for i = 1:NumExp
        tic;
        [wp, stats] = phipmactual(TF(k), Problem.A, ones(N, p+1),  err(k), sym,10);
        cpu1 = toc;
        cputimes(i) = cpu1;
    end
    tkchoicestate10(k,:) = stats;
    phipm10time(k) = mean(cputimes);
    phipm10err(k) = norm((wexact-wp)./wexact);
    
    for i = 1:NumExp
        tic;
        [wp, stats] = umphipm(TF(k), Problem.A, ones(N, p+1),err(k), sym);
        cpu1 = toc;
        cputimes(i) = cpu1;
    end
    stateum(k,:) = stats;
    umphipmtime(k) = mean(cputimes);
    umphipmerr(k) = norm((wexact-wp)./wexact);
    
    
    for i = 1:NumExp
        tic;
        [wp, stats] = uphipm1(TF(k), Problem.A, ones(N, p+1),err(k), sym);
        cpu1 = toc;
        cputimes(i) = cpu1;
    end
    tkchoicestate1(k,:) = stats;
    uphipm1time(k) = mean(cputimes);
    uphipm1err(k) = norm((wexact-wp)./wexact);
    
    for i = 1:NumExp
        tic;
        [wp, stats] = uphipm(TF(k), Problem.A, ones(N, p+1),err(k), sym);
        cpu1 = toc;
        cputimes(i) = cpu1;
    end
    tkchoicestateu(k,:) = stats;
    uphipmtime(k) = mean(cputimes);
    uphipmerr(k) = norm((wexact-wp)./wexact); 
end
%%
statephip
stateuphip
tkchoicestate
tkchoicestate10
stateum
tkchoicestate1
tkchoicestateu
timerecord = [phipmtime;phipm10time]

err = [phiperr ;uphiperr ;phipmerr ;umphipmerr ;uphipm1err;uphipmerr;phipm10err];

time = [phiptime./uphiptime ;phiptime./phipmtime; phiptime./umphipmtime ;phiptime./uphipm1time;phiptime./uphipmtime]


